
# library imports
library(tidyverse)
library(scales)
library(limma)
library(edgeR)
library(psych)

Warning message:
“replacing previous import ‘lifecycle::last_warnings’ by ‘rlang::last_warnings’ when loading ‘tibble’”
── Attaching packages ─────────────────────────────────────────────────────────────────────────────────────────── tidyverse 1.3.0 ──

✔ ggplot2 3.3.2     ✔ purrr   0.3.4
✔ tibble  3.0.3     ✔ dplyr   1.0.2
✔ tidyr   1.1.1     ✔ stringr 1.4.0
✔ readr   1.3.1     ✔ forcats 0.5.0

── Conflicts ────────────────────────────────────────────────────────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag()    masks stats::lag()


Attaching package: ‘scales’


The following object is masked from ‘package:purrr’:

    discard


The following object is masked from ‘package:readr’:

    col_factor


Attaching package: ‘psych’


The following objects are masked from ‘package:scales’:

    alpha, rescale


The following objects are masked from ‘package:ggplot2’:

    %+%, alpha


# ================== TMM normalization from DGEList object =====================
apply_tmm_factors <- function(y, color = NULL, plot = TRUE, table = TRUE, title = "TMM Normalized data") {
    # computes the tmm normalized data from the DGEList object
        # y - DGEList object
        # color - color vector for boxplots
        # plot - flag to show boxplot or not
        # table - flag to print norm factor tables or not
        # title - main title for boxplots
        # returns a dataframe with normalized intensities
    
    # compute and print "Sample loading" normalization factors
    lib_facs <- mean(y$samples$lib.size) / y$samples$lib.size
    if(table == TRUE) {
        cat("\nLibrary size factors:\n", 
            sprintf("%-5s -> %f\n", colnames(y$counts), lib_facs))
    }
    
    # compute and print TMM normalization factors
    tmm_facs <- 1/y$samples$norm.factors
    if(table == TRUE) {
        cat("\nTrimmed mean of M-values (TMM) factors:\n", 
            sprintf("%-5s -> %f\n", colnames(y$counts), tmm_facs))
    }
    
    # compute and print the final correction factors
    norm_facs <- lib_facs * tmm_facs
    if(table == TRUE) {
        cat("\nCombined (lib size and TMM) normalization factors:\n", 
            sprintf("%-5s -> %f\n", colnames(y$counts), norm_facs))
    }

    # compute the normalized data as a new data frame
    tmt_tmm <- as.data.frame(sweep(y$counts, 2, norm_facs, FUN = "*"))
    colnames(tmt_tmm) <- str_c(colnames(y$counts), "_tmm")
    
    # visualize results and return data frame
    if(plot == TRUE) {
        boxplot(log10(tmt_tmm), col = color, notch = TRUE, main = title)
    }
    tmt_tmm
}

# ============ runs an edgeR exact pairwise test ========================
pairwise_test <- function(delist, pair, p.value) {
    # runs an edgeR exact test and returns results
    et <- exactTest(delist, pair = pair)
    tt <- topTags(et, n = Inf, sort.by = "none")
    
    # print top tags
    cat("\n")
    print(topTags(et)$table, digits = 4)
    cat("\n")
    
    # count candidates
    print(summary(decideTestsDGE(et, p.value = p.value)))

    # see candidates on an MA plot
    plotMD(et, p.value = p.value)
    abline(h = c(-1, 1), col = "black")
    
    # check the p-value distribution
    pval <- ggplot(tt$table, aes(PValue)) + 
      geom_histogram(bins = 100, fill = "white", color = "black") + 
      geom_hline(yintercept = mean(hist(et$table$PValue, breaks = 100, 
                                    plot = FALSE)$counts[26:100])) +
      ggtitle("p-value distribution")
    print(pval) # this makes the plot show up
    
    tt # return top-tags object
}

# ================= reformat edgeR test results ================================
collect_results <- function(df, tt, x, xlab, y, ylab) {
    # Computes new columns and extracts some columns to make results frame
        # df - data in data.frame
        # tt - top tags table from edgeR test
        # x - columns for first condition
        # xlab - label for x
        # y - columns for second condition
        # ylab - label for y
        # returns a new dataframe
    
    # condition average vectors
    ave_x <- rowMeans(df[x])
    ave_y <- rowMeans(df[y])
    
    # FC, direction, candidates
    fc <- ifelse(ave_y > ave_x, (ave_y / ave_x), (-1 * ave_x / ave_y))
    direction <- ifelse(ave_y > ave_x, "up", "down")
    candidate <- cut(tt$FDR, breaks = c(-Inf, 0.01, 0.05, 0.10, 1.0), 
                     labels = c("high", "med", "low", "no"))
    
    # make data frame
    temp <- cbind(df[c(x, y)], data.frame(logFC = tt$logFC, FC = fc, 
                                          PValue = tt$PValue, FDR = tt$FDR, 
                                          ave_x = ave_x, ave_y = ave_y, 
                                          direction = direction, candidate = candidate, 
                                          Acc = tt$genes)) 
    
    # fix column headers for averages
    names(temp)[names(temp) %in% c("ave_x", "ave_y")]  <- str_c("ave_", c(xlab, ylab))    
    
    temp # return the data frame
}

# =============== log2_fold_changes ================================================
log2_fold_changes <- function(results, left, right, title) {
    # Makes log2 fold change plots by candidate
        # results - results data frame
        # left - lower FC log2 limit
        # right - upper FC log2 limit
        # title - plot title
    
    # see how many candidates by category
    cat("\n")
    print(results %>% count(candidate))
    cat("\n")

    # plot log2 fold-changes by category
    fc <- ggplot(results, aes(x = logFC, fill = candidate)) +
      geom_histogram(binwidth=0.1, color = "black") +
      facet_wrap(~candidate) +
      coord_cartesian(xlim = c(-left, right)) +
      ggtitle(title)
    print(fc)
}

# ========== Setup for MA and volcano plots ====================================
transform <- function(results, x, y) {
    # Make data frame with some transformed columns
        # results - results data frame
        # x - columns for x condition
        # y - columns for y condition
        # return new data frame
    df <- data.frame(log10((results[x] + results[y])/2), 
                     log2(results[y] / results[x]), 
                     results$candidate,
                     -log10(results$FDR))
    colnames(df) <- c("A", "M", "candidate", "P")
    
    df # return the data frame
}

# ========== MA plots using ggplot =============================================
MA_plots <- function(results, x, y, title) {
    # makes MA-plot DE candidate ggplots
        # results - data frame with edgeR results and some condition average columns
        # x - string for x-axis column
        # y - string for y-axis column
        # title - title string to use in plots
        # returns a list of plots 
    
    # uses transformed data
    temp <- transform(results, x, y)
    
    # 2-fold change lines
    ma_lines <- list(geom_hline(yintercept = 0.0, color = "black"),
                     geom_hline(yintercept = 1.0, color = "black", linetype = "dotted"),
                     geom_hline(yintercept = -1.0, color = "black", linetype = "dotted"))

    # make main MA plot
    ma <- ggplot(temp, aes(x = A, y = M)) +
        geom_point(aes(color = candidate, shape = candidate)) +
        scale_y_continuous(paste0("logFC (", y, "/", x, ")")) +
        scale_x_continuous("Ave_intensity") +
        ggtitle(title) + 
        ma_lines
    
    # make separate MA plots
    ma_facet <- ggplot(temp, aes(x = A, y = M)) +
        geom_point(aes(color = candidate, shape = candidate)) +
        scale_y_continuous(paste0("log2 FC (", y, "/", x, ")")) +
        scale_x_continuous("log10 Ave_intensity") +
        ma_lines +
        facet_wrap(~ candidate) +
        ggtitle(str_c(title, " (separated)"))

    # make the plots visible
    print(ma)
    print(ma_facet)
}    

# ========== Scatter plots using ggplot ========================================
scatter_plots <- function(results, x, y, title) {
    # makes scatter-plot DE candidate ggplots
        # results - data frame with edgeR results and some condition average columns
        # x - string for x-axis column
        # y - string for y-axis column
        # title - title string to use in plots
        # returns a list of plots
    
    # 2-fold change lines
    scatter_lines <- list(geom_abline(intercept = 0.0, slope = 1.0, color = "black"),
                          geom_abline(intercept = 0.301, slope = 1.0, color = "black", linetype = "dotted"),
                          geom_abline(intercept = -0.301, slope = 1.0, color = "black", linetype = "dotted"),
                          scale_y_log10(),
                          scale_x_log10())

    # make main scatter plot
    scatter <- ggplot(results, aes_string(x, y)) +
        geom_point(aes(color = candidate, shape = candidate)) +
        ggtitle(title) + 
        scatter_lines

    # make separate scatter plots
    scatter_facet <- ggplot(results, aes_string(x, y)) +
        geom_point(aes(color = candidate, shape = candidate)) +
        scatter_lines +
        facet_wrap(~ candidate) +
        ggtitle(str_c(title, " (separated)")) 

    # make the plots visible
    print(scatter)
    print(scatter_facet)
}

# ========== Volcano plots using ggplot ========================================
volcano_plot <- function(results, x, y, title) {
    # makes a volcano plot
        # results - a data frame with edgeR results
        # x - string for the x-axis column
        # y - string for y-axis column
        # title - plot title string
    
    # uses transformed data
    temp <- transform(results, x, y)
    
    # build the plot
    ggplot(temp, aes(x = M, y = P)) +
        geom_point(aes(color = candidate, shape = candidate)) +
        xlab("log2 FC") +
        ylab("-log10 FDR") +
        ggtitle(str_c(title, " Volcano Plot"))
}

# ============== individual protein expression plots ===========================
# function to extract the identifier part of the accesssion
get_identifier <- function(accession) {
    identifier <- str_split(accession, "\\|", simplify = TRUE)
    identifier[,3]
#    identifier <- accession
}

set_plot_dimensions <- function(width_choice, height_choice) {
    options(repr.plot.width=width_choice, repr.plot.height=height_choice)
}

plot_top_tags <- function(results, nleft, nright, top_tags, color = c("red", "blue")) {
    # results should have data first, then test results (two condition summary table)
    # nleft, nright are number of data points in each condition
    # top_tags is number of up and number of down top DE candidates to plot
    # get top ipregulated
    up <- results %>% 
        filter(logFC >= 0) %>%
        arrange(FDR)
    up <- up[1:top_tags, ]
    
    # get top down regulated
    down <- results %>% 
        filter(logFC < 0) %>%
        arrange(FDR)
    down <- down[1:top_tags, ]
    
    # pack them
    proteins <- rbind(up, down)
        
    color_bars = c(rep(color[1], nleft), rep(color[2], nright))
    for (row_num in 1:nrow(proteins)) {
        row <- proteins[row_num, ]
        vec <- as.vector(unlist(row[1:(nleft + nright)]))
        names(vec) <- colnames(row[1:(nleft + nright)])
        title <- str_c(get_identifier(row$Acc), ", int: ", scientific(mean(vec), 2), 
                       ", FDR: ", scientific(row$FDR, digits = 3), 
                       ", FC: ", round(row$FC, digits = 1),
                       ", ", row$candidate)
        barplot(vec, col = color_bars, main = title,
                cex.main = 1.0, cex.names = 0.7, cex.lab = 0.7)
    }    
}

# ============== CV function ===================================================
CV <- function(df) {
    # Computes CVs of data frame rows
        # df - data frame, 
        # returns vector of CVs (%)
    ave <- rowMeans(df)    # compute averages
    sd <- apply(df, 1, sd) # compute standard deviations
    cv <- 100 * sd / ave   # compute CVs in percent (last thing gets returned)
}


# load the IRS-normalized data and check the table
data_import <- read_tsv("labeled_grouped_protein_summary_TMT_9_IRS_normalized.txt", 
                        guess_max = 5854)

# "Filter" column flags contams and decoys
# "Missing" column flags proteins without reporter ion intensities (full sets missing)
# the table from pandas is sorted so the rows we want come first
data_all <- filter(data_import, is.na(Filter), is.na(Missing))
data_sl <- select(data_all, contains("SLNorm_"))
data_irs <- select(data_all, contains("IRSNorm_"))

# save gene names for edgeR so we can double check that results line up
accessions <- data_all$Accession

# see how many rows of data we have
length(accessions)

# set some larger plot dimensions
set_plot_dimensions(9, 9)

Parsed with column specification:
cols(
  .default = col_double(),
  Accession = col_character(),
  Identical = col_character(),
  Similar = col_character(),
  OtherLoci = col_character(),
  Filter = col_character(),
  Missing = col_character(),
  Description = col_character()
)

See spec(...) for full column specifications.


# define the groups (after IRS)
E15_irs <- select(data_irs, contains("_E15_"))
E18_irs <- select(data_irs, contains("_E18_"))
P0_irs <- select(data_irs, contains("_P0_"))
P3_irs <- select(data_irs, contains("_P3_"))
P6_irs <- select(data_irs, contains("_P6_"))
P9_irs <- select(data_irs, contains("_P9_"))

# put groups together into a single data frame
tmt_data <- bind_cols(E15_irs, E18_irs, P0_irs, P3_irs, P6_irs, P9_irs)

# Set color vector:
colors <- c(rep("blue", 3), rep("red", 6), rep("dark green", 6), rep("grey", 3))

# make indexes for young and old
E15 <- 1:3
mid <- 4:9
older <- 10:15
P9 <- 16:18


# get the biological sample data into a DGEList object
group <- c(rep('E15', 3), rep('mid', 6), rep('older', 6), rep('P9', 3))
y <- DGEList(counts = tmt_data, group = group, genes = accessions)

# run TMM normalization so we have the results in y
y <- calcNormFactors(y)

tmt_tmm <- apply_tmm_factors(y, table = TRUE, color = colors)

# check the clustering
plotMDS(y, col = colors, main = "Data after TMM")

Library size factors:
 IRSNorm_E15_1 -> 1.031476
 IRSNorm_E15_2 -> 1.008094
 IRSNorm_E15_3 -> 0.986400
 IRSNorm_E18_1 -> 1.010526
 IRSNorm_E18_2 -> 1.012694
 IRSNorm_E18_3 -> 0.997802
 IRSNorm_P0_1 -> 0.998726
 IRSNorm_P0_2 -> 0.991580
 IRSNorm_P0_3 -> 0.999595
 IRSNorm_P3_1 -> 0.986234
 IRSNorm_P3_2 -> 0.999611
 IRSNorm_P3_3 -> 1.005907
 IRSNorm_P6_1 -> 0.980859
 IRSNorm_P6_2 -> 0.993474
 IRSNorm_P6_3 -> 1.011724
 IRSNorm_P9_1 -> 0.974078
 IRSNorm_P9_2 -> 1.004781
 IRSNorm_P9_3 -> 1.009580

Trimmed mean of M-values (TMM) factors:
 IRSNorm_E15_1 -> 0.712297
 IRSNorm_E15_2 -> 0.660870
 IRSNorm_E15_3 -> 0.666675
 IRSNorm_E18_1 -> 0.839775
 IRSNorm_E18_2 -> 0.704975
 IRSNorm_E18_3 -> 0.849973
 IRSNorm_P0_1 -> 0.908618
 IRSNorm_P0_2 -> 0.896906
 IRSNorm_P0_3 -> 0.880678
 IRSNorm_P3_1 -> 1.077557
 IRSNorm_P3_2 -> 1.213998
 IRSNorm_P3_3 -> 1.069485
 IRSNorm_P6_1 -> 1.224558
 IRSNorm_P6_2 -> 1.314099
 IRSNorm_P6_3 -> 1.317240
 IRSNorm_P9_1 -> 1.463104
 IRSNorm_P9_2 -> 1.354746
 IRSNorm_P9_3 -> 1.501006

Combined (lib size and TMM) normalization factors:
 IRSNorm_E15_1 -> 0.734718
 IRSNorm_E15_2 -> 0.666219
 IRSNorm_E15_3 -> 0.657608
 IRSNorm_E18_1 -> 0.848614
 IRSNorm_E18_2 -> 0.713924
 IRSNorm_E18_3 -> 0.848105
 IRSNorm_P0_1 -> 0.907461
 IRSNorm_P0_2 -> 0.889354
 IRSNorm_P0_3 -> 0.880321
 IRSNorm_P3_1 -> 1.062723
 IRSNorm_P3_2 -> 1.213527
 IRSNorm_P3_3 -> 1.075802
 IRSNorm_P6_1 -> 1.201120
 IRSNorm_P6_2 -> 1.305523
 IRSNorm_P6_3 -> 1.332683
 IRSNorm_P9_1 -> 1.425177
 IRSNorm_P9_2 -> 1.361223
 IRSNorm_P9_3 -> 1.515386


# compute dispersions and plot BCV
y <- estimateDisp(y)
plotBCV(y, main = "BCV plot of IRS/TMM normalized data")

Design matrix not provided. Switch to the classic mode.


# use a 0.05 cutoff
tt = pairwise_test(y, c("E15", "mid"), 0.05)

                      genes   logFC logCPM    PValue       FDR
410   sp|P07356|ANXA2_MOUSE  1.2510  8.413 1.485e-19 6.311e-16
848    sp|P62317|SMD2_MOUSE -0.9485  6.796 3.067e-16 4.779e-13
402   sp|Q8BFX1|RN187_MOUSE  2.6494  8.140 3.806e-16 4.779e-13
1188  sp|O35658|C1QBP_MOUSE -0.8486  5.942 4.499e-16 4.779e-13
1374  sp|A2AWP0|BIRC7_MOUSE  1.9059  5.903 1.032e-14 8.771e-12
389    sp|Q64448|CXA3_MOUSE  1.3453  8.097 3.037e-14 2.151e-11
1376  sp|P05132|KAPCA_MOUSE  1.2097  5.591 5.065e-14 3.074e-11
615  tr|D3Z2H9|D3Z2H9_MOUSE -0.7327  7.462 1.989e-13 1.056e-10
166    sp|Q9DCL9|PUR6_MOUSE  0.9877  9.805 4.478e-13 2.114e-10
467   sp|Q91WJ8|FUBP1_MOUSE -0.8942  7.805 1.024e-12 4.352e-10

       mid-E15
Down       492
NotSig    3213
Up         544


# get the results summary
results <- collect_results(tmt_tmm, tt$table, E15, "E15", mid, "mid")

# make column names unique by adding comparison (for the accumulated frame)
results_temp  <- results
colnames(results_temp) <- str_c(colnames(results), "_E15_mid")

# accumulate testing results
all_results <- cbind(accessions, results_temp)


# count candidates and check log2 fold change distributions
log2_fold_changes(results, 3, 3, "E15 vs mid-age lenses logFC distributions by candidate")

  candidate    n
1      high  596
2       med  440
3       low  310
4        no 2903


# make MA plots
MA_plots(results, "ave_E15", "ave_mid", "E15 vs mid-age lenses")


# make scatter plots
scatter_plots(results, "ave_E15", "ave_mid", "E15 vs mid-age lenses")


# make a volcano plot
volcano_plot(results, "ave_E15", "ave_mid", "E15 vs mid-age lenses")


# look at the top candidates
set_plot_dimensions(9, 3.5)

n_left <- 3
n_right <- 6
top_tags <- 100
color <- c("blue", "red") 

plot_top_tags(results, n_left, n_right, top_tags, color)
set_plot_dimensions(9, 9)


# use a 0.05 cutoff
tt = pairwise_test(y, c("mid", "older"), 0.05)

                      genes   logFC logCPM    PValue       FDR
234  tr|E9Q616|E9Q616_MOUSE  1.3464  9.289 7.877e-35 3.347e-31
772   sp|Q8K0S2|ERIC5_MOUSE  0.9720  7.402 4.493e-18 9.546e-15
172   sp|Q62318|TIF1B_MOUSE -0.8207  9.319 7.985e-17 8.965e-14
366    sp|Q99MD9|NASP_MOUSE -1.1305  8.174 8.439e-17 8.965e-14
811   sp|Q3ULJ0|GPD1L_MOUSE  0.9674  6.807 8.122e-16 6.902e-13
1057  sp|Q920Q8|NS1BP_MOUSE  0.7774  6.075 2.120e-15 1.501e-12
1191  sp|Q9QYH6|MAGD1_MOUSE  0.7820  6.081 5.392e-14 3.273e-11
6     sp|P23927|CRYAB_MOUSE  1.9323 14.472 1.616e-13 7.715e-11
266   sp|P10107|ANXA1_MOUSE  1.0138  9.252 1.634e-13 7.715e-11
166    sp|Q9DCL9|PUR6_MOUSE  0.7786  9.805 3.349e-13 1.423e-10

       older-mid
Down         603
NotSig      3107
Up           539


# get the results summary
results <- collect_results(tmt_tmm, tt$table, mid, "mid", older, "older")

# make column names unique by adding comparison (for the accumulated frame)
results_temp  <- results
colnames(results_temp) <- str_c(colnames(results), "_mid_older")

# accumulate testing results
all_results <- cbind(accessions, results_temp)


# count candidates and check log2 fold change distributions
log2_fold_changes(results, 3, 3, "Mid-age vs Older-age lenses logFC distributions by candidate")

  candidate    n
1      high  718
2       med  424
3       low  285
4        no 2822


# make MA plots
MA_plots(results, "ave_mid", "ave_older", "Mid-age vs older-age lenses")


# make scatter plots
scatter_plots(results, "ave_mid", "ave_older", "Mid-age vs older-age lenses")


# make a volcano plot
volcano_plot(results, "ave_mid", "ave_older", "Mid-age vs older-age lenses")


# look at the top candidates
set_plot_dimensions(9, 3.5)

n_left <- 6
n_right <- 6
top_tags <- 100
color <- c("red", "dark green") 

plot_top_tags(results, n_left, n_right, top_tags, color)
set_plot_dimensions(9, 9)


write.table(all_results, "Khan-2018_results_3-age-groups_edgeR-exact_TMM.txt", sep = "\t",
           row.names = FALSE, na =  " ")


sessionInfo()

R version 3.5.3 (2019-03-11)
Platform: x86_64-apple-darwin15.6.0 (64-bit)
Running under: macOS  10.16

Matrix products: default
BLAS: /Library/Frameworks/R.framework/Versions/3.5/Resources/lib/libRblas.0.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/3.5/Resources/lib/libRlapack.dylib

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] psych_2.0.7     edgeR_3.24.3    limma_3.38.3    scales_1.1.1   
 [5] forcats_0.5.0   stringr_1.4.0   dplyr_1.0.2     purrr_0.3.4    
 [9] readr_1.3.1     tidyr_1.1.1     tibble_3.0.3    ggplot2_3.3.2  
[13] tidyverse_1.3.0

loaded via a namespace (and not attached):
 [1] Rcpp_1.0.5       locfit_1.5-9.4   lubridate_1.7.9  lattice_0.20-41 
 [5] assertthat_0.2.1 digest_0.6.25    IRdisplay_0.7.0  R6_2.4.1        
 [9] cellranger_1.1.0 repr_1.1.0       backports_1.1.8  reprex_0.3.0    
[13] evaluate_0.15    httr_1.4.2       pillar_1.4.6     rlang_1.0.2     
[17] uuid_0.1-4       readxl_1.3.1     rstudioapi_0.11  blob_1.2.1      
[21] labeling_0.3     splines_3.5.3    munsell_0.5.0    broom_0.7.0     
[25] compiler_3.5.3   modelr_0.1.8     pkgconfig_2.0.3  base64enc_0.1-3 
[29] mnormt_1.5-6     htmltools_0.4.0  tidyselect_1.1.0 crayon_1.3.4    
[33] dbplyr_1.4.4     withr_2.2.0      grid_3.5.3       nlme_3.1-148    
[37] jsonlite_1.7.0   gtable_0.3.0     lifecycle_0.2.0  DBI_1.1.0       
[41] magrittr_1.5     cli_3.3.0        stringi_1.4.6    farver_2.0.3    
[45] fs_1.5.0         xml2_1.3.2       ellipsis_0.3.2   generics_0.0.2  
[49] vctrs_0.4.1      IRkernel_1.1.1   tools_3.5.3      glue_1.6.2      
[53] hms_0.5.3        parallel_3.5.3   colorspace_1.4-1 rvest_0.3.6     
[57] pbdZMQ_0.3-3     haven_2.3.1

Age	Description	Biological Replicates	Pooled Mice Lenses
E15	Embryonic day 15	3	23 embryos
E18	Embryonic day 18	3	10 embryos
P0	Newborn	3	8 pups
P3	3-day old	3	8 pups
P6	6-day old	3	8 pups
P9	9-day old	3	8 pups

Channel	Plex 1	Plex 2	Plex 3
126C	E15	E15	E15
127N	E18	E18	E18
127C	P0	P0	P0
128N	P3	P3	P3
128C	P6	P6	P6
129N	P9	P9	P9
129C	Mix?	Mix?	Mix?
130N	unused	unused	unused
130C	unused	unused	unused
131N	unused	unused	unused

Khan 2018 Mouse Lens Re-analysis¶

Lens changes with age groups notebook¶

E15 -vs E18+P0 and E18+P0 -vs- P3+P6 in pairwise comparisons¶

TMM normalization is used after IRS here¶

Prepared by Phil Wilmarth, PSR Core, OHSU¶

December 26, 2023¶

Overview¶

Experimental Groups:¶

Sample Key:¶

This notebook will compare lenses of different ages¶

References¶

Load the necessary R libraries¶

Define common functions for notebook use¶

Read the IRS-normalized TMT intensity data¶

Organize the data frame by biological group¶

Load data into edgeR data structures¶

Groups separate left-to-right¶

EdgeR statistical testing starts here¶

Compute the shared variance trend¶

E15 lenses (n=3) versus E18+P0 lenses (n=6)¶

We have around 500 DE candidates in each direction¶

Compute averages, flag candidates, and save the pair-wise test results¶

Count candidates and look at fold-change distributions¶

Main summary plots¶

MA plots¶

The non-candidates and candidates have a typical pattern¶

Scatter plots¶

Volcano plot¶

Check some individual proteins¶

E18 and P0 seem similar¶

Mid-age lenses (n=6) versus older lenses (n=6)¶

We have around 600 DE candidates in each direction¶

Compute averages, flag candidates, and save the pair-wise test results¶

Count candidates and look at fold-change distributions¶

Main summary plots¶

MA plots¶

The non-candidates and candidates have a typical pattern¶

Scatter plots¶

Volcano plot¶

Check some individual proteins¶

Top candidates look convincing¶

Summary¶

Save the all_results frame to TSV file¶

Log the session information¶

Overview ¶

Read the IRS-normalized TMT intensity data ¶

Load data into edgeR data structures ¶

Summary ¶

Save the `all_results` frame to TSV file¶